Optimal Design of Minimum-Power Stimuli for Spiking Neurons 
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In this article, we study optimal control problems of spiking neurons whose dynamics are described 
by a phase model. We design minimum-power current stimuli (controls) that lead to targeted spiking 
times of neurons, where the cases with unbounded and bounded control amplitude are considered. 
We show that theoretically the spiking period of a neuron, modeled by phase dynamics, can be 
arbitrarily altered by a smooth control. However, if the control amplitude is bounded, the range 
of possible spiking times is constrained and determined by the bound, and feasible spiking times 
are optimally achieved by piecewise continuous controls. We present analytic expressions of these 
minimum-power stimuli for spiking neurons and illustrate the optimal solutions with numerical 
simulations. 

PACS numbers: 02.30.Yy, 87.19. lr, 87.19.11 



I. INTRODUCTION 

Control of neurons and hence the nervous system by 
external current stimuli (controls) has received increased 
scientific attention in recent years for its wide range of 
applications from deep brain stimulation to oscillatory 
neurocomputers [1-3] . Conventionally, neuron oscillators 
are represented by phase-reduced models, which form a 
standard nonlinear system [4, 5]. Intensive studies us- 
ing phase models have been carried out, for example, on 
the investigation of the patterns of synchrony that re- 
sult from the type and architecture of coupling [6, 7] and 
on the response of large groups of oscillators to external 
stimuli [8, 9], where the inputs to the neuron systems 
were initially defined and the dynamics of neural popu- 
lations were analyzed in detail. 

Recently, control theoretic approaches have been em- 
ployed to design external stimuli that drive neurons to 
behave in a desired way. For example, a multilinear feed- 
back control technique has been used to control the in- 
dividual phase relation between coupled oscillators [10]; 
a nonlinear feedback approach has been employed to en- 
gineer complex dynamic structures and synthesize deli- 
cate synchronization features of nonlinear systems [11]; 
and our recent work has illustrated controllability of a 
network of neurons with different natural oscillation fre- 
quencies adopting tools from geometric control theory 
[12]. 

There has been an increase in the demand for con- 
trolling not only the collective behavior of a network of 
oscillators but also the behavior of each individual oscil- 
lator. It is feasible to change the spiking periods of oscil- 
lators or tune the individual phase relationship between 
coupled oscillators by the use of electric stimuli [10, 13]. 
Minimum-power stimuli that elicit spikes of a neuron at 
specified times close to the natural spiking time were an- 
alyzed [8]. Optimal waveforms for the entrainment of 



weakly forced oscillators that maximize the locking range 
have been calculated, where first and second harmonics 
were used to approximate the phase response curve [14] . 
These optimal controls were found mainly based on the 
calculus of variations, which restricts the optimal solu- 
tions to the class of smooth controls and the bound of 
the control amplitude was not taken into account. 

In this paper, we apply the Pontryagin's maximum 
principle [15, 16] to derive minimum-power controls that 
spike a neuron at desired time instants. We consider both 
cases when the available control amplitude is unbounded 
and bounded. The latter is of practical importance due 
to physical limitations of experimental equipment and 
the safety margin for neurons, e.g., the requirement of 
a mild brain stimulations in neurological treatments for 
Parkinson's disease and epilepsy. 

This paper is organized as follows. In Section II, we 
introduce the phase model for spiking neurons and formu- 
late the related optimal control problem. In Section III, 
we derive minimum-power controls associated with spec- 
ified spiking times in the absence and presence of con- 
trol amplitude constraints, in which various phase mod- 
els including sinusoidal PRC, SNIPER PRC, and theta 
neuron models are considered. In addition, we present 
examples and simulations to demonstrate the resulting 
optimal control strategies. 



II. OPTIMAL CONTROL OF SPIKING 
NEURONS 

A periodically spiking or firing neuron can be consid- 
ered as a periodic oscillator governed by the nonlinear 
dynamical equation of the form 



Z(0)I(t), 



(i) 
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where 9 is the phase of the oscillation, f(0) and Z(8) 
are real-valued functions giving the neuron's baseline dy- 
namics and its phase response, respectively, and I(t) is an 
external current stimulus [4]. The nonlinear dynamical 
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system described in (1) is referred to as the phase model 
for the neuron. The assumption that Z(9) vanishes only 
on isolated points and that / (9) > are made so that 
a full revolution of the phase is possible. By convention, 
neuron spikes occur when 9 = 2mr, where n e N. In 
the absence of any input I(t), the neuron spikes period- 
ically at its natural frequency, while the periodicity can 
be altered in a desired manner by an appropriate choice 
of J(t). 

In this article, we study optimal design of neural inputs 
that lead to the spiking of neurons at a specified time T 
after spiking at time t = 0. In particular, we find the 
stimulus that fires a neuron with minimum power, which 
is formulated as the following optimal control problem, 



mm 

i(t) 

s.t. 



/ I(tf dt 
Jo 



(2) 



9 = f(9) + Z(9)I(t), 
9(0) = 0, 9(T) = 2tt 
\I(t)\<M, Vi, 



where M > is the amplitude bound of the current stim- 
ulus J(i). Note that instantaneous or arbitrarily delayed 
spiking of a neuron is possible if I(t) is unbounded, i.e., 
M — oo; however, the range of feasible spiking periods 
of a neuron described as in (1) is restricted with a finite 
M. We consider both unbounded and bounded cases. 



where A is the Lagrange multiplier. The necessary op- 
timality conditions according to the Maximum Principle 
give 



dH 

A = -— = -Xz d I cos 9, 



(5) 



dH 
dl 



and 
I satisfies 



21 + \zd.suv9 = 0. Hence, the optimal current 



I = --Xz d sm9. 



(6) 



Substituting (6) into (3) and (5), the optimal control 
problem is then transformed to a boundary value prob- 
lem, which characterizes the optimal trajectories of 9(t) 
and \(t), 



UJ - 



ZjX . 2/i 

-4- sin 2 9, 



z 2 X 2 

A = sin 9 cos 9, 



(7) 

(8) 



with boundary conditions 9(0) = and 9(T) = 2n while 
A(0) and A(T) are unspecified. 

Additionally, since the Hamiltonian is not explicitly 
dependent on time, the optimal triple (A, 9, I) satisfies 
H(X,9,I) = c, V0 < t < T, where c is a constant. To- 
gether with (6), this yields 



III. MINIMUM-POWER STIMULUS FOR 
SPECIFIED FIRING TIME 

We consider the minimum-power optimal control prob- 
lem of spiking neurons as formulated in (2) for various 
phase models including sinusoidal PRC, SNIPER PRC, 
and theta neuron. 



Sinusoidal PRC Phase Model 



Consider the sinusoidal PRC model, 
9 = uj + Zds'm9 ■ I(t), 



(3) 



where ui is the natural oscillation frequency of the neuron 
and Zd is a model-dependent constant. The neuron de- 
scribed by this phase model spikes periodically with the 
period T = 2it/uj in the absence of any external input, 
i.e., I(t) = 0. 



1. Spiking Neurons with Unbounded Control 

The optimal current profile can be derived by Pontrya- 
gin's Maximum Principle [15]. Given the optimal control 
problem as in (2), we form the control Hamiltonian 



H = I Z 



X(u! + Zd sin9 ■ I), 



(4) 



sin 2 9X 2 +wX = c, 



(9) 



and the constant c = cjAo is obtained by the initial con- 
ditions 9(0) = and A(0) = Ao, which is undetermined. 
Then, the optimal multiplier can be found by solving the 
above quadratic equation (9), which gives 



A = 



2ui ± 2 Juj 2 - luXqz 2 sin 2 6 



z 2 sin 2 9 



and the optimal trajectory of 9 follows 



= Tv^ 2 — wAoz^ sin 9, 



(10) 



(11) 



by plugging A in (10) into (7). Integrating (11) by separa- 
tion of variables, we find the spiking time T with respect 
to the initial condition Ao, 



1 



\J lj 2 — wAoz^ sin 2 8 



d9, (12) 



where F denotes the elliptic integral. Note that we 
choose the positive sign in (11) since the negative velocity 
indicates the backward phase evolution. Therefore, given 
a desired spiking time T of the neuron, the initial value, 
Ao, corresponding to the optimal trajectory of the mul- 
tiplier can be found via the one-to-one relation in (12). 
Consequently, the optimal trajectories of 9 and A can be 
easily computed by evolving (7) and (8) forward in time. 
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Phase (6) 



FIG. 1. Extremals of sinusoidal PRC system with Zd = 1 and 
u = 1 



Plugging (10) into (6), we obtain the optimal feedback 
law 



r = 



-co + <\J oj 2 — u}\qz^ sin 2 
Zd sin 6 



(13) 



which drives the neuron from 6(0) — to 6(T) — 2ir with 
minimum power. 

The feasibility of spiking the neuron at a desired time 
T largely depends on the initial value of the multiplier, 
Ao- It is clear from (12) and (11) that a complete 2ir rev- 
olution is impossible when Ao > uij 'z 2 . This fact also can 
be seen from FIG. 1, where the system evolution defined 
by (7) and (8) for Zd = 1 and to = 1 with respect to dif- 
ferent Ao values (9 = axis) is illustrated. When Ao = 0, 
according to equation (12), the spiking period is equal to 
the natural spiking period, 2tt/uj, and no external stim- 
ulus needs to be applied, i.e., I*(t) = 0, Vt € [0, 2w/u]. 
Since, from (12), T is a monotonically increasing func- 
tion of Ao for fixed w and Zd and, from (11), the average 
phase velocity decreases when Ao increases, the spiking 
time T > 2-k/uj for A > and T < 2it/uj for A < 0. 
FIG. 2 shows variation of the spiking time T with the Ao 
corresponding to the optimal trajectories for different u 
values with 2^ = 1. 

The relation between the spiking time T and required 
minimum power E = min I 2 (t)dt is evident via a sim- 
ple sensitivity analysis [17]. Since a small change in the 
initial condition, d9, and a small change in the initial 
time, dt, result in a small change in power according to 
dE = X(t)d6 - H(t)d6, it follows that [17] 



dE 

"at 



= H = c = wAo 



(14) 



T < 2-njuj (T > 2-tt/lj), we see that the required mini- 
mum power increases if we move away from the natural 
spiking time. 

The minimum-power stimulus I* as in (13) plotted 
with respect to time and phase for various spiking times 
T = 3, 5, 10, 12 with uj = 1 and Zd = 1 are shown in 
FIG. 3(a) and 3(b), respectively. The respective optimal 
trajectories of \(6) and 0(t) for these spiking times are 
presented in FIG. 3(c) and 3(d). 



2. Spiking Neurons with Bounded Control 

In practice, the amplitude of stimuli in physical sys- 
tems are limited, so we consider spiking the sinusoidal 
neuron with bounded control amplitude, namely, in the 
optimal control problem (2), \I(t)\ < M < 00 for all 
t € [0, T], where T is the desired spiking period. In this 
case, there exists a range of feasible spiking periods de- 
pending on the value of M, in contrast to the previous 
case where any desired spiking time is feasible. We first 
observe that given this bound M, the minimum time it 
takes to spike a neuron can be achieved by choosing the 
control that keeps the phase velocity 8 maximum over 
t G [0,T]. Such a time-optimal control, for Zd > 0, can 
be characterized by a switching, i.e., 



Train 



M 
-M 



for 
for 



< 9 < 7T 
7T < 6 < 2tt 



(15) 



Consequently, the spiking time with I^, r 
puted using (3) and (15), which yields 



can be com- 



T 



M 



2tt - 4 tan" 1 \z d M/ ^-z'jM 2 + w 2 | 



v- 



(16) 



It follows that J*, derived in (13), is the minimum-power 
stimulus that spikes the neuron at a desired spiking time 
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This implies that E increases with initial time t for Ao < 
and decreases for Ao > 0. Since the increment of the ini- 
tial time is equivalent to the decrement of spiking time 
T, 8E/8T = ujXq. Since A < (Ao > 0) corresponds to 



FIG. 2. Variation of the spiking time, T, with respect to the 
initial multiplier value, Ao, leading to optimal trajectories, 
with different values of cj and z<j = 1 for sinusoidal PRC 
model. 
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(c) (d) 



FIG. 3. Optimal solutions for various spiking times T — 
3, 5, 10, 12 for sinusoidal PRC model with Zd = 1 and u = 1. 
(a) The minimum-power control I* . (b) Variation of I* with 
phase 6. (c) Variation of the optimal multiplier, A*, with 9. 
(d) Optimal phase trajectories following 7*. 

T if |7*| < M for all t £ [0,T]. However, there exists 
a shortest possible spiking time by I* given the bound 
M. Simple first and second order optimality conditions 
applied to (13) find that the maximum value of I* occurs 
at 8 = it/ 2 for A < and at 8 = 3tt/2 for A > 0. There- 
fore, the Ao for the shortest spiking time with control I* 
satisfying |/*(t)| < M can be calculated by substituting 
I* = M and 8 = ir/2 to the equation (13), and then from 
(12) we obtain this shortest spiking period 

f 2rr 1 
T 1 = / -^^==^^^^^=^^= (17) 

rmn / I \ J 

J ° yjuj 2 + z d M(z d M + 2u) sin 2 (9) 

Note that T^ in < T^ m . According to (3) when M > 
uj/zd, arbitrarily large spiking times can be achieved by 
making 9 arbitrary close to zero. Therefore we consider 
two cases for M > oj/zj, and M < uj/zd- 

Case I: M > uj/zd- Since I* takes the maximum 
value at 9 — Zir/2 for Ao > 0, we have |/*| < (to — 
y/cu 2 — ujAqZj)/ Zd, which leads to |7*| < ujjzd < M for 
Ao > 0. This implies that I* is the minimum-power 
control for any desired spiking time T > 2tt/uj when 
M > uj/zd, and hence for any spiking time T > Tj^^. 
Variation of the maximum value of the control I* with 
spiking time T for uj — 1 and Zd = 1 is depicted in FIG. 4. 
Shorter spiking times T £ ^Tmim^min) are feasible but, 
due to the bound M, can not be achieved by I* since 
it requires a control with amplitude greater than M for 
some t £ [0, T]. However, these spiking times can be op- 
timally achieved by applying controls switching between 
I* and I* min . 

Let the desired spiking time T £ [T^ in ,T^ in ). Then, 



there exist two angles 8\ = sin 1 [— 2Muj / (zdM 2 + 
ZdUjXo)] and 92 = it — 9\ where I* meets the bound M. 
When 9 £ (9 U 9 2 ), I* > M and we take 1(9) = M for 9 £ 
[6*i, 6*2]. The Hamiltonian of the system when 9 £ [6\, 6*2] 
is then, from (4), H = M 2 +X(oj+z d sin 9 M). If the triple 
(A, 9, M) is optimal, then if is a constant, which gives 
A = (H - M 2 )/(u + z d M sin 9). This multiplier satisfies 
the adjoint equation (5), and therefore 1(9) = M is op- 
timal for 9 £ [81,82]- Similarly, by symmetry, I* < —AI 
when 9 £ [83, 6*4] , where #3 = n + 9\ and 84 — 2ir — 9\ , if 
the desired spiking time T £ [T^ in , T^ in ). It can be eas- 
ily shown by the same fashion that 1(9) — —M is optimal 
in the interval 8 £ [83, 6*4] . 

Therefore, the minimum-power optimal control that 
spikes the neuron at T £ [T^ in , T^ in ) can be character- 
ized by four switchings between I* and M, i.e., 

{r < 9 < 0i 

M 9 1 <9<9 2 

I* 9 2 <9<9 3 (18) 

—M 9 3 <9<9 4: 

I* 9 4 <9< 2tt. 

The initial value of the multiplier, Ao, resulting in the 
optimal trajectory, can then be found according to the 
desired spiking time T £ [T^ in ,T^ in ) through the rela- 
tion 

T= / = d8+ / —dd. 

J yjui 2 -uj\ a zlsin 2 6 J 6l UJ + ZdM sin (8) 

FIG. 5(a) shows the relation between Ao and T by If 
for M = 2.5, Zd — 1, and ui = 1. From (16) the minimum 
possible spiking time with this control bound M = 2.5 is 
Tmi n = 2.735 and from (17) the minimum spiking time 
by I* is T^ in = 3.056. Thus, in this example, any desired 
spiking time T > 3.056 can be optimally achieved by I* 
whereas any T £ [2.735, 3.056) can be optimally obtained 
by If as in (18). FIG. 5(b) illustrates the bounded and 
unbounded optimal controls that fire the neuron at T = 



5 




2.5 2.8 



FIG. 5. (a) Variation of the spiking time T £ PnSm^min) 
for sinusoidal PRC model with initial multiplier value, Ao, 
for the bound of control amplitude M = 2.5. (b) Minimum- 
power controls with (M = 2.5) and without a constraint on 
the control amplitude for sinusoidal PRC model with T — 2.8, 
Zd = 1, and uj = 1. 
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(a) (b) 

FIG. 6. (a) Variation of the spiking time T £ (T»,oii^«ii] 
with the initial value of the multiplier, Ao, for sinusoidal PRC 
model when M = 0.55. (b) Minimum-power controls with 
(M — 0.55) and without a constraint on the control amplitude 
for sinusoidal PRC model with T — 10, za = 1, and u = 1. 



2.8, where I* is the minimum-power stimulus when the 
control amplitude is not limited and I* is the minimum- 
power stimulus when the bound M = 2.5. I* drives the 
neuron from 9(0) = to 9(2.8) = 2vr with 13.54 units of 
power whereas /* requires 14.13 units. 

Case II : M < uo/z d . In contrast with Case I in the 
previous section, achieving arbitrarily large spiking times 
is not feasible with a bound M < ujj z d . In this case, the 
longest possible spiking time is achieved by 



Traax 



-M for < 9 < 7T, 
M for 7T < 9 < 2tt. 



The spiking time of the neuron under this control is, 



rpM _ 
max 



27T + 4 tan 



z d M/^-z 2 M 2 +uj< 



y/-z* d M 2 + uj 2 



(19) 



and the longest spiking time feasible with control I* is 
given by 



T 1 



^Juj 2 + z d M(z d M - 2uj) sin 2 (6») 



(20) 



Then, by similar analysis for Case I, any spiking time 
T *= Km^mm) for a given M < uj/z d can be achieved 
with the minimum-power control /* as given in (18), any 
T G [T^in' ^max] can achieved with minimum power 
by I* in (13), and moreover any T <E (T^ ax , T^ ax ] can be 
obtained by switching between I* and I* max - The corre- 
sponding switching angles are 65 = sin _1 [2A/o;/(z ( ;M 2 + 
z d uj\o)}, #6 = 7T — 05, 67 = 7T+#5 and 6^ = 2-k — 6$, and the 
minimum-power optimal control for T € (T/ 
characterized by 



< 9 < 9 5 
9 5 <9<9 6 
9 6 <9 <9 7 
7 <9<9$ 
9 8 <9< 2vr. 



T M 1 is 

max 5 maxi ^ 




The Ao resulting in the optimal trajectory by /| can be 
calculated according to the given T G (T^ r '' 
the relation 



T 1 via 

max'' max! 



T : 



\j uJ 1 — luXqZj sin 2 $ 



=d9 



q k uj — ZdM sin( 



-<m. 



FIG. 6(a) shows the relation between Ao and T by 
I* for M = 0.55, z d = 1, and uj = 1. From (19) 
the maximum possible spiking time with M — 0.55 is 
Tmax — 10.312 and from (20) the maximum spiking time 
feasible by I* is I max = 9.006. Therefore, in this ex- 
ample, any desired spiking time T G (9.006, 10.312] can 
be obtained with minimum power by the use of I^ ■ FIG. 
6(b) illustrates the bounded and unbounded optimal con- 
trols that spike the neuron at T = 10, where I* is the 
minimum-power stimulus when the control amplitude is 
not limited and is the minimum-power stimulus when 
M = 0.55. I* drives the neuron from 9(0) = to 
9(10) = 2ir with 2.193 units of power whereas J| requires 
2.327 units. 

A summary of the optimal (minimum-power) spiking 
scenarios for a prescribed spiking time of the neuron gov- 
erned by the sinusoidal phase model (3) is illustrated in 
FIG. 7. 



B. SNIPER PRC and Theta Neuron Phase Models 

We now consider the SNIPER PRC model in which 
f(9) = uj and Z(9) = z d (l — cos9), where z d > and 
uj > 0. That is, 



uj + z d (l - cos 9)I(t). 



(21) 



The minimum-power stimuli for spiking neurons modeled 
by this phase model can be easily derived with analogous 
analysis described previously in III A 1 and III A 2 for the 
sinusoidal PRC phase model. 
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FIG. 7. A summery of the optimal control strategies for the 
sinusoidal PRC neuron for (a) M > uj/z d , (b) M < uj/z d . 
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FIG. 8. Optimal solutions for various spiking times T — 
3,5,10,12 for Sniper PRC model with z d — 1 and cj = 1. 
(a) The minimum-power control J*, (b) Variation of J* with 
phase 9. (c) Variation of the optimal multiplier, A*, with 6. 
(d) Optimal phase trajectories following I*. 



1. Spiking Neurons with Unbounded Control 

Employing the maximum principle as in III A 1 , the 
minimum-power stimulus that spikes the SNIPER neuron 
at a desired time T can be derived and given by 



I* 



—lj + \/uj 2 — (jJ\qZj(1 — cost 
z d (l — cos 9) 



(22) 



where Ao corresponding to the optimal trajectory is de- 
termined through the integral relation with T, 



T = 



271 



1 



o \Jw 2 — wAqz 2 (1 — cos( 



-A9. 



The minimum-power stimuli I* plotted with respect to 
time and phase for various spiking times T = 3, 5, 10, 12 



with parameter values z d =l and u> = 1 are illustrated 
in FIG. 8(a) and 8(b), respectively. The corresponding 
optimal trajectories of X(6) and 9{t) for these spiking 
times are displayed in FIG. 8(c) and 8(d). 



2. Spiking Neurons with Bounded Control 

When the amplitude of the available stimulus is lim- 
ited, i.e., \I(t)\ < M, the control that achieves the 
shortest spiking time for the SNIPER neuron modeled 
in (21) is given by I% min = M > for < 6 < 2tt, 
since 1 — cos^ > for all 9 € [0, 2tt],. As a result, 
the shortest possible spiking time with this control is 
T" m = 2tt/Vw 2 + 2z d ujM. Also, the shortest spiking 
time achieved by the control I* in (22) given the bound 
M is given by 



T r = 

min 



l 



v/w 2 + z d M(z d M + lo)(1 



(23) 



Similar to the sinusoidal PRC case, the longest possible 
spiking time of the neuron varies with the control bound 
M. If M > u>/(2zd), an arbitrarily large spiking time 
is achievable, however, if M < uj/(2zd) there exists a 
maximum spiking time. 

Case I: M >u>/{2z d ). Any spiking time T £ [T£ in , oo) 
is possible with control /* but a shorter spiking time T € 
[ T min' T min) requires switching between I* and 
which is characterized by two switchings, 



n 



0<9 <6i 
9 1 <9<2tt-9 1 
2tt - 6 X < 9 < 2?r 



(24) 



where 6± = cos -1 [l + 2ujM/(z d M 2 + z d uj\ )] . The ini- 
tial value Ao which results in the optimal trajectory is 
given by, 



T = 



\/ to 2 — ljXqz^ cos G) 2 



--de 



1 ut-\-Zfi A'/(l — cos 6 



-de. 



FIG. 9(a) illustrates the relation between Ao and T £ 
[Tl n , T% in ) by It for M = 2,z d = 1, and u = 1. In this 
case, the shortest feasible spiking time is = 2.09 

and the shortest with the control /* is T^ in — 3.18. Any 
spiking time in the interval (2.09, 3.18] is achievable by Jj* 
in (24) with minimum-power. FIG. 9(b) illustrates the 
unbounded and bounded, with M = 2, optimal stimuli 
that fire the neuron at T — 3 with minimum-power. 

Case II: M < uj/(2z d ). In this case there exists a 
longest possible spiking time which is achieved by I^ax = 
—M for all 6 G [0, 2tt]. The longest spiking time feasible 
with the control I* as in (22) is given by 



T 



I* 



271 



1 



o v 7 ^ + z d M(z d M -2oj)(1 



Therefore, any spiking time T € [T min ,T min ) for a given 
M < lj I z d can be achieved with the minimum-power 
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FIG. 9. (a) Variation of the spiking time T G [T„, 
with the initial multiplier value, Ao, for SNIPER PRC model 
when M = 2. (b) Minimum-power controls with (M = 2) and 
without a constraint on the control amplitude for SNIPER 
PRC model with T = 3, z d = 1, and w = 1. 



T £j can be 



control 7* as given in (24), any T G [T^ in , 
achieved with minimum power by I* in (22), and more- 
over any T G (T^^, T^ ax ] can be obtained by switching 



2 tan -1 U/lb tan ((^> — 71") /2)] +7r, the theta neuron model 
can be transformed to a SNIPER PRC with z d = lj/2 
in (21). Therefore, all of the results developed for the 
SNIPER PRC can be directly applied to the theta neu- 
ron phase model. 
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FIG. 10. (a) Variation of the spiking time T G (Tm„,r m ni] 
with the initial multiplier value, Ao, for SNIPER PRC model 
when M = 0.3. (b) Minimum-power controls with (M = 
0.3) and with out a constraint on the control amplitude for 
SNIPER PRC model with T = 9.8, z d = 1, and w = 1. 



between I* and It 



that is, 




o < e < 2 

&2 < < 2tt - 6» 2 
2tt - 6 2 < 9 < 2tt 



where 2 = cos" 1 [l - 2ojM/(z d M 2 + z d oj\ )]. The A 
associated with the optimal trajectory is determined via 
the relation with the desired spiking time T, 



T = 



JO 



yj w 2 -cjA q z^(1-c 



-.do 



uj — ZdM(l—c 



-de. 



FIG. 10(a) illustrates the relation between Ao and T G 
(TL^Zxl by II for M = 0.3, z d = 1, and u = 1. 
In this case, the longest feasible spiking time is T^ ax = 
9.935 and the longest with the control I* is T^ ax = 8.596. 
The unbounded and bounded, with M = 0.3, optimal 
stimuli that fire the neuron at T = 9.8 with minimum- 
power are illustrated in FIG. 10(b). 

A summary of the optimal (minimum-power) spik- 
ing scenarios for a prescribed spiking time of the neu- 
ron governed by the SNIPER PRC model in (21) can 
be illustrated analogously to FIG. 7(a) and 7(b) for 
M > uj/(2zd) and M < w/(2zd), respectively. 

The theta neuron phase model is described by the dy- 
namical equation 

9 = 1 + cos6> + z d (l - cos0)(/(i) (25) 

where If, > is the baseline current and the nat- 
ural frequency ui of the theta neuron is given by 
2\flh. In fact, by a coordinate transformation, 



IV. CONCLUSION AND FUTURE WORK 

In this paper, we studied various phase-reduced models 
that describe the dynamics of a neuron system. We con- 
sidered the design of minimum-power stimuli for spiking a 
neuron at a specified time instant and formulated this as 
an optimal control problem. We investigated both cases 
when the control amplitude is unbounded and bounded, 
for which we found analytic expressions of optimal feed- 
back control laws. In particular for the bounded control 
case, we characterized the ranges of possible spiking peri- 
ods in terms of the control bound. Moreover, minimum- 
power stimuli for steering any nonlinear oscillators of the 
form as in (1) between desired states can be derived fol- 
lowing the steps presented in this article. In addition, 
the charge-balanced constraint [18] can be readily incor- 
porated into this framework as well. 

The optimal control of a single neuron system inves- 
tigated in this work illustrates the fundamental limit of 
spiking a neuron with external stimuli and provides a 
benchmark structure that enables us to study optimal 
control of a spiking neural network with different individ- 
ual oscillation frequencies. Our recent work [12] proved 
that simultaneous spiking of a network of neurons is pos- 
sible; however, optimal control of such a spiking neu- 
ral network has not been studied. We finally note that 
although one-dimensional phase models are reasonably 
accurate to describe the dynamics of neurons, study- 
ing higher dimensional models such as that of Hodgkin- 
Huxley is essential for more accurate computation of op- 
timal neural inputs. 
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